
# ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
# DESCRIPTION
# ______________________________________________________________________________

# This code creates Figures 4, 6, B.1, C.1–6, and D.1.

# Preamble ---------------------------------------------------------------------

rm(list = ls())

# Packages ---------------------------------------------------------------------

library(tidyverse)
library(ggpmisc)
library(ggpubr)
library(grid)
library(haven)
library(patchwork)
library(signs)

# Main Analysis (ZIOP) ---------------------------------------------------------

GeomSegment$draw_key <-  function (data, params, linewidth) {
  draw_key_vpath <- function (data, params, linewidth) {
    segmentsGrob(0.1,
                 0.5,
                 0.9,
                 0.5,
                 gp = gpar(
                   col = alpha(data$colour, data$alpha),
                   lwd = data$linewidth * .pt,
                   lty = data$linetype
                 ))
  }
  grobTree(draw_key_vpath(data, params, linewidth),
           draw_key_point(data, params, linewidth))
}

## Data ------------------------------------------------------------------------

main_rebel <- read_dta("./results/main_rebel.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "rebelabductever" ~ "Abduction",
      parm == "drugs" ~ 'Contraband',
      parm == "terrcontrol" ~ "Territorial Control",
      parm == "leftist" ~ "Leftist Ideology",
      parm == "frontlineprevbest" ~ "Women Combatants",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "statesv_max" ~ "State SV"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

main_state <- read_dta("./results/main_state.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "statepressgangever" ~ "Press-Ganging",
      parm == "lntroopqual" ~ "ln(Troop Quality)",
      parm == "fcstate" ~ "Women Combatants",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "rebelsv_max" ~ "Rebel SV"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

main_rebel$new_parm <- 
  factor(main_rebel$new_parm, levels = rev(unique(main_rebel$new_parm)))

main_state$new_parm <- 
  factor(main_state$new_parm, levels = rev(unique(main_state$new_parm)))

## Plot ------------------------------------------------------------------------

space_between_bars <- 0.1

main_rebel$type_y_adj =
  scale(as.numeric(as.factor(main_rebel$model_name))) * space_between_bars

main_rebel$y =
  as.numeric(as.numeric(as.factor(main_rebel$new_parm)) + main_rebel$type_y_adj)

main_rebel_plot <- main_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  annotate(
    "rect",
    xmin = -0.45,
    xmax = 0.45,
    ymin = 12.55,
    ymax = 13.45,
    linewidth = 0.2,
    fill = NA,
    color = "black"
  ) +
  annotate(
    "segment",
    x = 0.45,
    xend = 2,
    y = 13,
    yend = 13,
    linewidth = 0.2
  ) +
  annotate(
    "rect",
    xmin = -0.45,
    xmax = 0.45,
    ymin = 2.55,
    ymax = 3.45,
    linewidth = 0.2,
    fill = NA,
    color = "black"
  ) +
  annotate(
    "segment",
    x = 0.45,
    xend = 2,
    y = 3,
    yend = 3,
    linewidth = 0.2
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  facet_wrap(paste0("Rebel SV\nObservations = 626") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:15,
    labels = rev(unique(main_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    legend.position = "none",
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off")

main_rebel_plot_inset1 <- main_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 3,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:15,
    labels = rev(unique(main_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    plot.background = element_rect(fill = NA, color = NA),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(xlim = c(-0.3, 0.3), ylim = c(13, 13))

main_rebel_plot_inset2 <- main_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 3,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:15,
    labels = rev(unique(main_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    plot.background = element_rect(fill = NA, color = NA),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(xlim = c(-0.05, 0.05), ylim = c(3, 3))

space_between_bars <- 0.1

main_state$type_y_adj =
  scale(as.numeric(as.factor(main_state$model_name))) * space_between_bars

main_state$y =
  as.numeric(as.numeric(as.factor(main_state$new_parm)) + main_state$type_y_adj)

main_state_plot <- main_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  annotate(
    "rect",
    xmin = -0.48,
    xmax = 0.48,
    ymin = 2.55,
    ymax = 3.45,
    linewidth = 0.2,
    fill = NA,
    color = "black"
  ) +
  annotate(
    "segment",
    x = -0.48,
    xend = -2,
    y = 3,
    yend = 3,
    linewidth = 0.2
  ) +
  annotate(
    "rect",
    xmin = -0.48,
    xmax = 0.48,
    ymin = 1.55,
    ymax = 2.45,
    linewidth = 0.2,
    fill = NA,
    color = "black"
  ) +
  annotate(
    "segment",
    x = 0.48,
    xend = 2,
    y = 2,
    yend = 2,
    linewidth = 0.2
  ) +
  guides(
    color = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    shape = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), byrow = TRUE)
  ) +
  facet_wrap(paste0("State SV\nObservations = 620") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:13,
    labels = rev(unique(main_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    legend.direction = "vertical",
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(12.5, 3, 0, 3, "pt"),
    legend.spacing.x = unit(15, "pt"),
    legend.spacing.y = unit(-5, "pt"),
    legend.key.width = unit(25, "pt"),
    legend.key.size = unit(12.5, "pt"),
    legend.title = element_text(size = 10, hjust = 0.5),
    legend.text = element_text(size = 9),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off")

main_state_plot_inset_1 <- main_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 3,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:13,
    labels = rev(unique(main_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.direction = "vertical",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    legend.title = element_blank(),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key.width = unit(20, "pt"),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2),
    legend.position = "none",
    plot.background = element_rect(fill = NA, color = NA)
  ) +
  coord_cartesian(xlim = c(-0.02, 0.02), ylim = c(3, 3))

main_state_plot_inset_2 <- main_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 3,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:13,
    labels = rev(unique(main_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.direction = "vertical",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    legend.title = element_blank(),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key.width = unit(20, "pt"),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2),
    legend.position = "none",
    plot.background = element_rect(fill = NA, color = NA)
  ) +
  coord_cartesian(xlim = c(-0.1, 0.1), ylim = c(2, 2))

main_rebel_plot2 <- main_rebel_plot +
  annotation_custom(
    ggplotGrob(main_rebel_plot_inset1),
    xmin = 0.8,
    xmax = 3.7,
    ymin = 11.725,
    ymax = 13.725
  ) +
  annotation_custom(
    ggplotGrob(main_rebel_plot_inset2),
    xmin = 0.8,
    xmax = 3.7,
    ymin = 1.7,
    ymax = 3.7
  )

main_state_plot2 <- main_state_plot +
  annotation_custom(
    ggplotGrob(main_state_plot_inset_1),
    xmin = -0.95,
    xmax = -4.15,
    ymin = 1.775,
    ymax = 3.675
  ) +
  annotation_custom(
    ggplotGrob(main_state_plot_inset_2),
    xmin = 0.85,
    xmax = 4.05,
    ymin = 0.775,
    ymax = 2.675
  )

leg <- get_legend(main_state_plot)

leg_plot <- as_ggplot(leg)

design <- "
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  13
  13
"

main_plot <- main_rebel_plot2 + main_state_plot2 + leg_plot +
  plot_layout(design = design) &
  theme(legend.position = "none",
        plot.margin = margin(3, 3, 1.75, 2.5, "pt"))

ggsave(
  plot = main_plot,
  "figures/coefficient plots/main.pdf",
  device = cairo_pdf,
  width = 10,
  height = 6.5
)

# Main Analysis (OP and ZIOP) --------------------------------------------------

## Data ------------------------------------------------------------------------

main_rebel_op <- read_dta("./results/main_rebel_op.dta") %>%
  filter(!parm %in% c("cut1", "cut2", "cut3")) %>%
  mutate(model_name = case_when(eq == "rebelsv_max" ~ "OP"))

main_rebel_ziop <- read_dta("./results/main_rebel.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(model_name = case_when(
    eq == "rebelsv_max" ~ "ZIOP (Outcome)",
    eq == "inflate" ~ "ZIOP (Inflation)"
  ))

main_rebel <- rbind(main_rebel_op, main_rebel_ziop)

main_rebel <- main_rebel %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "rebelabductever" ~ "Abduction",
      parm == "drugs" ~ 'Contraband',
      parm == "terrcontrol" ~ "Territorial Control",
      parm == "leftist" ~ "Leftist Ideology",
      parm == "frontlineprevbest" ~ "Women Combatants",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "statesv_max" ~ "State SV"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

main_rebel$new_parm <- 
  factor(main_rebel$new_parm, levels = rev(unique(main_rebel$new_parm)))

main_rebel$model_name <- 
  factor(main_rebel$model_name, levels = rev(unique(main_rebel$model_name)))

main_state_op <- read_dta("./results/main_state_op.dta") %>%
  filter(!parm %in% c("cut1", "cut2", "cut3")) %>%
  mutate(model_name = case_when(eq == "statesv_max" ~ "OP"))

main_state_ziop <- read_dta("./results/main_state.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(model_name = case_when(
    eq == "statesv_max" ~ "ZIOP (Outcome)",
    eq == "inflate" ~ "ZIOP (Inflation)"
  ))

main_state <- rbind(main_state_op, main_state_ziop)

main_state <- main_state %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "statepressgangever" ~ "Press-Ganging",
      parm == "lntroopqual" ~ "ln(Troop Quality)",
      parm == "fcstate" ~ "Women Combatants",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "rebelsv_max" ~ "Rebel SV"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

main_state$new_parm <- 
  factor(main_state$new_parm, levels = rev(unique(main_state$new_parm)))

main_state$model_name <- 
  factor(main_state$model_name, levels = rev(unique(main_state$model_name)))

## Plot ------------------------------------------------------------------------

space_between_bars <- 0.1

main_rebel$type_y_adj =
  scale(as.numeric(as.factor(main_rebel$model_name))) * space_between_bars

main_rebel$y =
  as.numeric(as.numeric(as.factor(main_rebel$new_parm)) + main_rebel$type_y_adj)

main_rebel_plot <- main_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  annotate(
    "rect",
    xmin = -0.45,
    xmax = 0.45,
    ymin = 12.55,
    ymax = 13.45,
    linewidth = 0.2,
    fill = NA,
    color = "black"
  ) +
  annotate(
    "segment",
    x = 0.45,
    xend = 2,
    y = 13,
    yend = 13,
    linewidth = 0.2
  ) +
  annotate(
    "rect",
    xmin = -0.45,
    xmax = 0.45,
    ymin = 2.55,
    ymax = 3.45,
    linewidth = 0.2,
    fill = NA,
    color = "black"
  ) +
  annotate(
    "segment",
    x = 0.45,
    xend = 2,
    y = 3,
    yend = 3,
    linewidth = 0.2
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), order = 2)
  ) +
  facet_wrap(paste0("Rebel SV\nObservations = 626") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:15,
    labels = rev(unique(main_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22, 21)) +
  scale_color_manual(values = c("blue", "red", "limegreen")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    legend.position = "none",
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off")

main_rebel_plot_inset1 <- main_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 3,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:15,
    labels = rev(unique(main_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22, 21)) +
  scale_color_manual(values = c("blue", "red", "limegreen")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    plot.background = element_rect(fill = NA, color = NA),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(xlim = c(-0.3, 0.3), ylim = c(13, 13))

main_rebel_plot_inset2 <- main_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 3,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:15,
    labels = rev(unique(main_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22, 21)) +
  scale_color_manual(values = c("blue", "red", "limegreen")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    plot.background = element_rect(fill = NA, color = NA),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(xlim = c(-0.05, 0.05), ylim = c(3, 3))

space_between_bars <- 0.1

main_state$type_y_adj =
  scale(as.numeric(as.factor(main_state$model_name))) * space_between_bars

main_state$y =
  as.numeric(as.numeric(as.factor(main_state$new_parm)) + main_state$type_y_adj)

main_state_plot <- main_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  annotate(
    "rect",
    xmin = -0.48,
    xmax = 0.48,
    ymin = 2.55,
    ymax = 3.45,
    linewidth = 0.2,
    fill = NA,
    color = "black"
  ) +
  annotate(
    "segment",
    x = -0.48,
    xend = -2,
    y = 3,
    yend = 3,
    linewidth = 0.2
  ) +
  annotate(
    "rect",
    xmin = -0.48,
    xmax = 0.48,
    ymin = 1.55,
    ymax = 2.45,
    linewidth = 0.2,
    fill = NA,
    color = "black"
  ) +
  annotate(
    "segment",
    x = 0.48,
    xend = 2,
    y = 2,
    yend = 2,
    linewidth = 0.2
  ) +
  guides(
    color = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    shape = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), byrow = TRUE)
  ) +
  facet_wrap(paste0("State SV\nObservations = 620") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:13,
    labels = rev(unique(main_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22, 21)) +
  scale_color_manual(values = c("blue", "red", "limegreen")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    legend.direction = "vertical",
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(15, 3, 0, 3, "pt"),
    legend.spacing.x = unit(15, "pt"),
    legend.spacing.y = unit(-5, "pt"),
    legend.key.width = unit(25, "pt"),
    legend.key.size = unit(11, "pt"),
    legend.title = element_text(size = 10, hjust = 0.5),
    legend.text = element_text(size = 9),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off")

main_state_plot

main_state_plot_inset_1 <- main_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 3,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:13,
    labels = rev(unique(main_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22, 21)) +
  scale_color_manual(values = c("blue", "red", "limegreen")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.direction = "vertical",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    legend.title = element_blank(),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key.width = unit(20, "pt"),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2),
    legend.position = "none",
    plot.background = element_rect(fill = NA, color = NA)
  ) +
  coord_cartesian(xlim = c(-0.02, 0.02), ylim = c(3, 3))

main_state_plot_inset_2 <- main_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 3,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:13,
    labels = rev(unique(main_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22, 21)) +
  scale_color_manual(values = c("blue", "red", "limegreen")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.direction = "vertical",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    legend.title = element_blank(),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key.width = unit(20, "pt"),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2),
    legend.position = "none",
    plot.background = element_rect(fill = NA, color = NA)
  ) +
  coord_cartesian(xlim = c(-0.1, 0.1), ylim = c(2, 2))

main_rebel_plot2 <- main_rebel_plot +
  annotation_custom(
    ggplotGrob(main_rebel_plot_inset1),
    xmin = 0.8,
    xmax = 3.7,
    ymin = 11.725,
    ymax = 13.725
  ) +
  annotation_custom(
    ggplotGrob(main_rebel_plot_inset2),
    xmin = 0.8,
    xmax = 3.7,
    ymin = 1.7,
    ymax = 3.7
  )

main_state_plot2 <- main_state_plot +
  annotation_custom(
    ggplotGrob(main_state_plot_inset_1),
    xmin = -0.95,
    xmax = -4.15,
    ymin = 1.775,
    ymax = 3.675
  ) +
  annotation_custom(
    ggplotGrob(main_state_plot_inset_2),
    xmin = 0.85,
    xmax = 4.05,
    ymin = 0.775,
    ymax = 2.675
  )

leg <- get_legend(main_state_plot)

leg_plot <- as_ggplot(leg)

design <- "
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  13
  13
"

main_plot <- main_rebel_plot2 + main_state_plot2 + leg_plot +
  plot_layout(design = design) &
  theme(legend.position = "none",
        plot.margin = margin(3, 3, 1.75, 2.5, "pt"))

ggsave(
  plot = main_plot,
  "figures/coefficient plots/main_op.pdf",
  device = cairo_pdf,
  width = 10,
  height = 6.5
)

# Robustness Check 1 -----------------------------------------------------------

## Data ------------------------------------------------------------------------

rc1_rebel <- read_dta("./results/rc1_rebel.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "rebelabductever" ~ "Abduction",
      parm == "drugs" ~ 'Contraband',
      parm == "terrcontrol" ~ "Territorial Control",
      parm == "leftist" ~ "Leftist Ideology",
      parm == "frontlineprevbest" ~ "Women Combatants",
      parm == "exlgender" ~ "Exclusion by Gender",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "statesv_max" ~ "State SV"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

rc1_state <- read_dta("./results/rc1_state.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "statepressgangever" ~ "Press-Ganging",
      parm == "lntroopqual" ~ "ln(Troop Quality)",
      parm == "fcstate" ~ "Women Combatants",
      parm == "exlgender" ~ "Exclusion by Gender",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "rebelsv_max" ~ "Rebel SV"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

rc1_rebel$new_parm <- 
  factor(rc1_rebel$new_parm, levels = rev(unique(rc1_rebel$new_parm)))

rc1_state$new_parm <- 
  factor(rc1_state$new_parm, levels = rev(unique(rc1_state$new_parm)))

## Plot ------------------------------------------------------------------------

space_between_bars <- 0.1

rc1_rebel$type_y_adj =
  scale(as.numeric(as.factor(rc1_rebel$model_name))) * space_between_bars

rc1_rebel$y =
  as.numeric(as.numeric(as.factor(rc1_rebel$new_parm)) + rc1_rebel$type_y_adj)

rc1_rebel_plot <- rc1_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  facet_wrap(paste0("Rebel SV\nObservations = 626") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:15,
    labels = rev(unique(rc1_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) + coord_cartesian(clip = "off")

space_between_bars <- 0.1

rc1_state$type_y_adj =
  scale(as.numeric(as.factor(rc1_state$model_name))) * space_between_bars

rc1_state$y =
  as.numeric(as.numeric(as.factor(rc1_state$new_parm)) + rc1_state$type_y_adj)

rc1_state_plot <- rc1_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    shape = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), byrow = TRUE)
  ) +
  facet_wrap(paste0("State SV\nObservations = 620") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:13,
    labels = rev(unique(rc1_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    legend.direction = "vertical",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(12.5, 3, 0, 3, "pt"),
    legend.spacing.x = unit(15, "pt"),
    legend.spacing.y = unit(-5, "pt"),
    legend.key.width = unit(25, "pt"),
    legend.key.size = unit(12.5, "pt"),
    legend.title = element_text(size = 10, hjust = 0.5),
    legend.text = element_text(size = 9),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) + coord_cartesian(clip = "off")

leg <- get_legend(rc1_state_plot)

leg_plot <- as_ggplot(leg)

design <- "
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  13
  13
"

rc1_plot <- rc1_rebel_plot + rc1_state_plot + leg_plot +
  plot_layout(design = design) &
  theme(legend.position = "none",
        plot.margin = margin(3, 3, 1.75, 2.5, "pt"))

ggsave(
  plot = rc1_plot,
  "figures/coefficient plots/rc1.pdf",
  device = cairo_pdf,
  width = 10,
  height = 6.5
)

# Robustness Check 2 -----------------------------------------------------------

## Data ------------------------------------------------------------------------

rc2_rebel <- read_dta("./results/rc2_rebel.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3", "yrsp2", "yrsp3")) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "rebelabductever" ~ "Abduction",
      parm == "drugs" ~ 'Contraband',
      parm == "terrcontrol" ~ "Territorial Control",
      parm == "leftist" ~ "Leftist Ideology",
      parm == "frontlineprevbest" ~ "Women Combatants",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "statesv_max" ~ "State SV"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  ) %>%
  add_row(new_parm = "Year Splines",
          sig = 1,
          model_name = "ZIOP (Outcome)")

rc2_state <- read_dta("./results/rc2_state.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3"),
         str_detect(parm, "year") != TRUE) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "statepressgangever" ~ "Press-Ganging",
      parm == "lntroopqual" ~ "ln(Troop Quality)",
      parm == "fcstate" ~ "Women Combatants ",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "rebelsv_max" ~ "Rebel SV"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  ) %>%
  add_row(new_parm = "Year Dummies",
          sig = 1,
          model_name = "ZIOP (Outcome)")

rc2_rebel$new_parm <- 
  factor(rc2_rebel$new_parm, levels = rev(unique(rc2_rebel$new_parm)))

rc2_state$new_parm <- 
  factor(rc2_state$new_parm, levels = rev(unique(rc2_state$new_parm)))

## Plot ------------------------------------------------------------------------

space_between_bars <- 0.1

rc2_rebel$type_y_adj =
  scale(as.numeric(as.factor(rc2_rebel$model_name))) * space_between_bars

rc2_rebel$y =
  as.numeric(as.numeric(as.factor(rc2_rebel$new_parm)) + rc2_rebel$type_y_adj)

rc2_rebel_plot <- rc2_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  annotate(
    "text",
    x = 0,
    y = 1,
    label = paste("🗸"),
    color = "black",
    alpha = 1,
    size = 8
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  facet_wrap(paste0("Rebel SV\nObservations = 626") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:15,
    labels = rev(unique(rc2_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) + coord_cartesian(clip = "off", ylim = c(1, 15))

space_between_bars <- 0.1

rc2_state$type_y_adj =
  scale(as.numeric(as.factor(rc2_state$model_name))) * space_between_bars

rc2_state$y =
  as.numeric(as.numeric(as.factor(rc2_state$new_parm)) + rc2_state$type_y_adj)

rc2_state_plot <- rc2_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  annotate(
    "text",
    x = 0,
    y = 1,
    label = paste("🗸"),
    color = "black",
    alpha = 1,
    size = 8
  ) +
  annotate(
    "rect",
    xmin = -1,
    xmax = 1,
    ymin = 12.55,
    ymax = 13.45,
    linewidth = 0.2,
    fill = NA,
    color = "black"
  ) +
  annotate(
    "segment",
    x = 1,
    xend = 3,
    y = 13,
    yend = 13,
    linewidth = 0.2
  ) +
  guides(
    color = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    shape = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), byrow = TRUE)
  ) +
  facet_wrap(paste0("State SV\nObservations = 620") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:13,
    labels = rev(unique(rc2_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    legend.direction = "vertical",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(12.5, 3, 0, 3, "pt"),
    legend.spacing.x = unit(15, "pt"),
    legend.spacing.y = unit(-5, "pt"),
    legend.key.width = unit(25, "pt"),
    legend.key.size = unit(12.5, "pt"),
    legend.title = element_text(size = 10, hjust = 0.5),
    legend.text = element_text(size = 9),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off", ylim = c(1, 13))

rc2_state_plot_inset <- rc2_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  annotate(
    "text",
    x = 0,
    y = 1,
    label = paste("🗸"),
    color = "black",
    alpha = 1,
    size = 8
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 3,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:13,
    labels = rev(unique(rc2_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    plot.background = element_rect(fill = NA, color = NA),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(xlim = c(-1, 1), ylim = c(13, 13))

rc2_state_plot2 <- rc2_state_plot +
  annotation_custom(
    ggplotGrob(rc2_state_plot_inset),
    xmin = 2,
    xmax = 10,
    ymin = 11.55,
    ymax = 13.45
  )

leg <- get_legend(rc2_state_plot)

leg_plot <- as_ggplot(leg)

design <- "
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  13
  13
"

rc2_plot <- rc2_rebel_plot + rc2_state_plot2 + leg_plot +
  plot_layout(design = design) &
  theme(legend.position = "none",
        plot.margin = margin(3, 3, 1.75, 2.5, "pt"))

ggsave(
  plot = rc2_plot,
  "figures/coefficient plots/rc2.svg",
  width = 10,
  height = 6.5
) # Save as .pdf using Inkscape: File > Save As... > PDF

# Robustness Check 3 -----------------------------------------------------------

## Data ------------------------------------------------------------------------

rc3_rebel <- read_dta("./results/rc3_rebel.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "rebelabductever" ~ "Abduction",
      parm == "drugs" ~ 'Contraband',
      parm == "terrcontrol" ~ "Territorial Control",
      parm == "leftist" ~ "Leftist Ideology",
      parm == "frontlineprevbest" ~ "Women Combatants",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "statesv_max" ~ "State SV",
      parm == "l_rebelsv_max" ~ "Rebel SV (t \u2013 1)"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

rc3_state <- read_dta("./results/rc3_state.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "statepressgangever" ~ "Press-Ganging",
      parm == "lntroopqual" ~ "ln(Troop Quality)",
      parm == "fcstate" ~ "Women Combatants",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "rebelsv_max" ~ "Rebel SV",
      parm == "l_statesv_max" ~ "State SV (t \u2013 1)"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

rc3_rebel$new_parm <- 
  factor(rc3_rebel$new_parm, levels = rev(unique(rc3_rebel$new_parm)))

rc3_state$new_parm <- 
  factor(rc3_state$new_parm, levels = rev(unique(rc3_state$new_parm)))

## Plot ------------------------------------------------------------------------

space_between_bars <- 0.1

rc3_rebel$type_y_adj =
  scale(as.numeric(as.factor(rc3_rebel$model_name))) * space_between_bars

rc3_rebel$y =
  as.numeric(as.numeric(as.factor(rc3_rebel$new_parm)) + rc3_rebel$type_y_adj)

rc3_rebel_plot <- rc3_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  facet_wrap(paste0("Rebel SV\nObservations = 558") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:16,
    labels = rev(unique(rc3_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) + coord_cartesian(clip = "off")

space_between_bars <- 0.1

rc3_state$type_y_adj =
  scale(as.numeric(as.factor(rc3_state$model_name))) * space_between_bars

rc3_state$y =
  as.numeric(as.numeric(as.factor(rc3_state$new_parm)) + rc3_state$type_y_adj)

rc3_state_plot <- rc3_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    shape = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), byrow = TRUE)
  ) +
  facet_wrap(paste0("State SV\nObservations = 542") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:14,
    labels = rev(unique(rc3_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    legend.direction = "vertical",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(12.5, 3, 0, 3, "pt"),
    legend.spacing.x = unit(15, "pt"),
    legend.spacing.y = unit(-5, "pt"),
    legend.key.width = unit(25, "pt"),
    legend.key.size = unit(12.5, "pt"),
    legend.title = element_text(size = 10, hjust = 0.5),
    legend.text = element_text(size = 9),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) + coord_cartesian(clip = "off")

leg <- get_legend(rc3_state_plot)

leg_plot <- as_ggplot(leg)

design <- "
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  13
  13
"

rc3_plot <- rc3_rebel_plot + rc3_state_plot + leg_plot +
  plot_layout(design = design) &
  theme(legend.position = "none",
        plot.margin = margin(3, 3, 1.75, 2.5, "pt"))

ggsave(
  plot = rc3_plot,
  "figures/coefficient plots/rc3.pdf",
  device = cairo_pdf,
  width = 10,
  height = 6.75
)

# Robustness Check 4 -----------------------------------------------------------

## Data ------------------------------------------------------------------------

rc4_rebel <- read_dta("./results/rc4_rebel.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3"),
         str_detect(parm, "0b") != TRUE) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "rebelabductever" ~ "Abduction",
      parm == "drugs" ~ 'Contraband',
      parm == "terrcontrol" ~ "Territorial Control",
      parm == "leftist" ~ "Leftist Ideology",
      parm == "1.frontlineprevbest" ~ "Women Combatants: Occasional",
      parm == "2.frontlineprevbest" ~ "Women Combatants: Low",
      parm == "3.frontlineprevbest" ~ "Women Combatants: Moderate",
      parm == "4.frontlineprevbest" ~ "Women Combatants: High",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "1.statesv_max" ~ "State SV: Level 1",
      parm == "2.statesv_max" ~ "State SV: Level 2",
      parm == "3.statesv_max" ~ "State SV: Level 3",
      
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

rc4_state <- read_dta("./results/rc4_state.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3"),
         str_detect(parm, "0b") != TRUE) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "statepressgangever" ~ "Press-Ganging",
      parm == "lntroopqual" ~ "ln(Troop Quality)",
      parm == "2.fcstate" ~ "Women Combatants: High",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "1.rebelsv_max" ~ "Rebel SV: Level 1",
      parm == "2.rebelsv_max" ~ "Rebel SV: Level 2",
      parm == "3.rebelsv_max" ~ "Rebel SV: Level 3"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

rc4_rebel$new_parm <- 
  factor(rc4_rebel$new_parm, levels = rev(unique(rc4_rebel$new_parm)))

rc4_state$new_parm <- 
  factor(rc4_state$new_parm, levels = rev(unique(rc4_state$new_parm)))

## Plot ------------------------------------------------------------------------

space_between_bars <- 0.1

rc4_rebel$type_y_adj =
  scale(as.numeric(as.factor(rc4_rebel$model_name))) * space_between_bars

rc4_rebel$y =
  as.numeric(as.numeric(as.factor(rc4_rebel$new_parm)) + rc4_rebel$type_y_adj)

rc4_rebel_plot <- rc4_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  facet_wrap(paste0("Rebel SV\nObservations = 626") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:19,
    labels = rev(unique(rc4_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) + coord_cartesian(clip = "off")

space_between_bars <- 0.1

rc4_state$type_y_adj =
  scale(as.numeric(as.factor(rc4_state$model_name))) * space_between_bars

rc4_state$y =
  as.numeric(as.numeric(as.factor(rc4_state$new_parm)) + rc4_state$type_y_adj)

rc4_state_plot <- rc4_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    shape = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), byrow = TRUE)
  ) +
  facet_wrap(paste0("State SV\nObservations = 620") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 5,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:14,
    labels = rev(unique(rc4_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    legend.direction = "vertical",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(12.5, 3, 0, 3, "pt"),
    legend.spacing.x = unit(15, "pt"),
    legend.spacing.y = unit(-5, "pt"),
    legend.key.width = unit(25, "pt"),
    legend.key.size = unit(12.5, "pt"),
    legend.title = element_text(size = 10, hjust = 0.5),
    legend.text = element_text(size = 9),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) + coord_cartesian(clip = "off")

leg <- get_legend(rc4_state_plot)

leg_plot <- as_ggplot(leg)

design <- "
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  13
  13
  13
  13
  13
  13
"

rc4_plot <- rc4_rebel_plot + rc4_state_plot + leg_plot +
  plot_layout(design = design) &
  theme(legend.position = "none",
        plot.margin = margin(3, 3, 1.75, 2.5, "pt"))

ggsave(
  plot = rc4_plot,
  "figures/coefficient plots/rc4.pdf",
  device = cairo_pdf,
  width = 10,
  height = 7
)

# Robustness Check 5 -----------------------------------------------------------

## Data ------------------------------------------------------------------------

rc5_rebel <- read_dta("./results/rc5_rebel.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "rebelabductever" ~ "Abduction",
      parm == "drugs" ~ 'Contraband',
      parm == "terrcontrol" ~ "Territorial Control",
      parm == "leftist" ~ "Leftist Ideology",
      parm == "frontlineprevbest" ~ "Women Combatants",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "statesv_max" ~ "State SV"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

rc5_state <- read_dta("./results/rc5_state.dta") %>%
  filter(!parm %in% c(
    "_cons",
    "cut1",
    "cut2",
    "cut3",
    "yrsp1",
    "yrsp2",
    "yrsp3",
    "yrsp4"
  )) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "statepressgangever" ~ "Press-Ganging",
      parm == "lntroopqual" ~ "ln(Troop Quality)",
      parm == "fcstate" ~ "Women Combatants ",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "rebelsv_max" ~ "Rebel SV"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

rc5_rebel$new_parm <- 
  factor(rc5_rebel$new_parm, levels = rev(unique(rc5_rebel$new_parm)))

rc5_state$new_parm <- 
  factor(rc5_state$new_parm, levels = rev(unique(rc5_state$new_parm)))

## Plot ------------------------------------------------------------------------

space_between_bars <- 0.1

rc5_rebel$type_y_adj =
  scale(as.numeric(as.factor(rc5_rebel$model_name))) * space_between_bars

rc5_rebel$y =
  as.numeric(as.numeric(as.factor(rc5_rebel$new_parm)) + rc5_rebel$type_y_adj)

rc5_rebel_plot <- rc5_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  facet_wrap(paste0("Rebel SV\nObservations = 761") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:11,
    labels = rev(unique(rc5_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) + coord_cartesian(clip = "off", ylim = c(1, 11))

space_between_bars <- 0.1

rc5_state$type_y_adj =
  scale(as.numeric(as.factor(rc5_state$model_name))) * space_between_bars

rc5_state$y =
  as.numeric(as.numeric(as.factor(rc5_state$new_parm)) + rc5_state$type_y_adj)

rc5_state_plot <- rc5_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  annotate(
    "rect",
    xmin = -0.75,
    xmax = 0.75,
    ymin = 8.55,
    ymax = 9.45,
    linewidth = 0.2,
    fill = NA,
    color = "black"
  ) +
  annotate(
    "segment",
    x = 0.75,
    xend = 4,
    y = 9,
    yend = 9,
    linewidth = 0.2
  ) +
  guides(
    color = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    shape = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), byrow = TRUE)
  ) +
  facet_wrap(paste0("State SV\nObservations = 697") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:9,
    labels = rev(unique(rc5_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    legend.direction = "vertical",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(12.5, 3, 0, 3, "pt"),
    legend.spacing.x = unit(15, "pt"),
    legend.spacing.y = unit(-5, "pt"),
    legend.key.width = unit(25, "pt"),
    legend.key.size = unit(12.5, "pt"),
    legend.title = element_text(size = 10, hjust = 0.5),
    legend.text = element_text(size = 9),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) + coord_cartesian(clip = "off", ylim = c(1, 9))

rc5_state_plot_inset <- rc5_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 3,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:9,
    labels = rev(unique(rc5_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    plot.background = element_rect(fill = NA, color = NA),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(xlim = c(-0.6, 0.6), ylim = c(9, 9))

rc5_state_plot2 <- rc5_state_plot +
  annotation_custom(
    ggplotGrob(rc5_state_plot_inset),
    xmin = 1.425,
    xmax = 8,
    ymin = 7.7,
    ymax = 9.45
  )

leg <- get_legend(rc5_state_plot)

leg_plot <- as_ggplot(leg)

design <- "
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  13
  13
"

rc5_plot <- rc5_rebel_plot + rc5_state_plot2 + leg_plot +
  plot_layout(design = design) &
  theme(legend.position = "none",
        plot.margin = margin(3, 3, 1.75, 2.5, "pt"))

ggsave(
  plot = rc5_plot,
  "figures/coefficient plots/rc5.pdf",
  device = cairo_pdf,
  width = 10,
  height = 5
)

# Robustness Check 6 -----------------------------------------------------------

## Data ------------------------------------------------------------------------

rc6_rebel <- read_dta("./results/rc6_rebel.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "rebelabductever" ~ "Abduction",
      parm == "drugs" ~ 'Contraband',
      parm == "leftist" ~ "Leftist Ideology",
      parm == "frontlineprevbest" ~ "Women Combatants",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "statesv_max" ~ "State SV"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

rc6_state <- read_dta("./results/rc6_state.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "statepressgangever" ~ "Press-Ganging",
      parm == "lntroopqual" ~ "ln(Troop Quality)",
      parm == "fcstate" ~ "Women Combatants",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "rebelsv_max" ~ "Rebel SV"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

rc6_rebel$new_parm <- 
  factor(rc6_rebel$new_parm, levels = rev(unique(rc6_rebel$new_parm)))

rc6_state$new_parm <- 
  factor(rc6_state$new_parm, levels = rev(unique(rc6_state$new_parm)))

## Plot ------------------------------------------------------------------------

space_between_bars <- 0.1

rc6_rebel$type_y_adj =
  scale(as.numeric(as.factor(rc6_rebel$model_name))) * space_between_bars

rc6_rebel$y =
  as.numeric(as.numeric(as.factor(rc6_rebel$new_parm)) + rc6_rebel$type_y_adj)

rc6_rebel_plot <- rc6_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  facet_wrap(paste0("Rebel SV\nObservations = 104") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:13,
    labels = rev(unique(rc6_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) + coord_cartesian(clip = "off")

space_between_bars <- 0.1

rc6_state$type_y_adj =
  scale(as.numeric(as.factor(rc6_state$model_name))) * space_between_bars

rc6_state$y =
  as.numeric(as.numeric(as.factor(rc6_state$new_parm)) + rc6_state$type_y_adj)

rc6_state_plot <- rc6_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    shape = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), byrow = TRUE)
  ) +
  facet_wrap(paste0("State SV\nObservations = 134") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:12,
    labels = rev(unique(rc6_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    legend.direction = "vertical",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(12.5, 36, 0, -29, "pt"),
    legend.spacing.x = unit(5, "pt"),
    legend.spacing.y = unit(-5, "pt"),
    legend.key.width = unit(25, "pt"),
    legend.key.size = unit(12.5, "pt"),
    legend.title = element_text(size = 9, hjust = 0.5),
    legend.text = element_text(size = 8),
    legend.title.position = "left",
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) + coord_cartesian(clip = "off")

leg <- get_legend(rc6_state_plot)

leg_plot <- as_ggplot(leg)

design <- "
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  13
"

rc6_plot <- rc6_rebel_plot + rc6_state_plot + leg_plot +
  plot_layout(design = design) &
  theme(legend.position = "none",
        plot.margin = margin(3, 3, 1.75, 2.5, "pt"))

ggsave(
  plot = rc6_plot,
  "figures/coefficient plots/rc6.pdf",
  device = cairo_pdf,
  width = 10,
  height = 6
)

# Longitudinal Analysis --------------------------------------------------------

## Data ------------------------------------------------------------------------

longitudinal_rebel <- read_dta("./results/longitudinal_rebel.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict_duration" ~ "Secessionist Conflict ×\nConflict Duration",
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "rebelabductever" ~ "Abduction",
      parm == "drugs" ~ 'Contraband',
      parm == "terrcontrol" ~ "Territorial Control",
      parm == "leftist" ~ "Leftist Ideology",
      parm == "frontlineprevbest" ~ "Women Combatants",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "statesv_max" ~ "State SV"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

longitudinal_state <- read_dta("./results/longitudinal_state.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict_duration" ~ "Secessionist Conflict ×\nConflict Duration",
      parm == "secessionistconflict" ~ "Secessionist Conflict",
      parm == "relrebelstr" ~ "Relative Rebel Strength",
      parm == "intervention" ~ "Intervention Balance",
      parm == "statepressgangever" ~ "Press-Ganging",
      parm == "lntroopqual" ~ "ln(Troop Quality)",
      parm == "fcstate" ~ "Women Combatants",
      parm == "gender" ~ "Women Political Empowerment",
      parm == "regimetype" ~ "Regime Type",
      parm == "lngdppc" ~ "ln(GDPpc)",
      parm == "lnpop" ~ "ln(Population)",
      parm == "conflictintst" ~ "Conflict Intensity",
      parm == "duration" ~ "Conflict Duration",
      parm == "rebelsv_max" ~ "Rebel SV"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  )

longitudinal_rebel$new_parm <- 
  factor(longitudinal_rebel$new_parm, levels = rev(unique(longitudinal_rebel$new_parm)))

longitudinal_state$new_parm <- 
  factor(longitudinal_state$new_parm, levels = rev(unique(longitudinal_state$new_parm)))

## Plot ------------------------------------------------------------------------

space_between_bars <- 0.1

longitudinal_rebel$type_y_adj =
  scale(as.numeric(as.factor(longitudinal_rebel$model_name))) * space_between_bars

longitudinal_rebel$y =
  as.numeric(as.numeric(as.factor(longitudinal_rebel$new_parm)) + longitudinal_rebel$type_y_adj)

longitudinal_rebel_plot <- longitudinal_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  annotate(
    "rect",
    xmin = -0.5,
    xmax = 0.5,
    ymin = 15.55,
    ymax = 16.45,
    linewidth = 0.2,
    fill = NA,
    color = "black"
  ) +
  annotate(
    "segment",
    x = -0.5,
    xend = -2,
    y = 16,
    yend = 16,
    linewidth = 0.2
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  facet_wrap(paste0("Rebel SV\nObservations = 626") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:16,
    labels = rev(unique(longitudinal_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off", ylim = c(1, 16.3))

longitudinal_rebel_plot_inset <- longitudinal_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 3,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:16,
    labels = rev(unique(longitudinal_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    plot.background = element_rect(fill = NA, color = NA),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(xlim = c(-0.2, 0.2), ylim = c(16, 16))

space_between_bars <- 0.1

longitudinal_state$type_y_adj =
  scale(as.numeric(as.factor(longitudinal_state$model_name))) * space_between_bars

longitudinal_state$y =
  as.numeric(as.numeric(as.factor(longitudinal_state$new_parm)) + longitudinal_state$type_y_adj)

longitudinal_state_plot <- longitudinal_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  annotate(
    "rect",
    xmin = -0.575,
    xmax = 0.575,
    ymin = 13.55,
    ymax = 14.45,
    linewidth = 0.2,
    fill = NA,
    color = "black"
  ) +
  annotate(
    "segment",
    x = 0.575,
    xend = 2,
    y = 14,
    yend = 14,
    linewidth = 0.2
  ) +
  guides(
    color = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    shape = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), byrow = TRUE)
  ) +
  facet_wrap(paste0("State SV\nObservations = 620") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:14,
    labels = rev(unique(longitudinal_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    legend.direction = "vertical",
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(12.5, 3, 0, 3, "pt"),
    legend.spacing.x = unit(15, "pt"),
    legend.spacing.y = unit(-5, "pt"),
    legend.key.width = unit(25, "pt"),
    legend.key.size = unit(12.5, "pt"),
    legend.title = element_text(size = 10, hjust = 0.5),
    legend.text = element_text(size = 9),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off", ylim = c(1, 14.3))

longitudinal_state_plot_inset <- longitudinal_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(
    limits = symmetric_limits,
    n.breaks = 3,
    labels = signs_format(
      format = function(x)
        as.character(x)
    )
  ) +
  scale_y_continuous(
    breaks = 1:14,
    labels = rev(unique(longitudinal_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    plot.background = element_rect(fill = NA, color = NA),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(xlim = c(-0.1, 0.1), ylim = c(14, 14))

longitudinal_rebel_plot2 <- longitudinal_rebel_plot +
  annotation_custom(
    ggplotGrob(longitudinal_rebel_plot_inset),
    xmin = -5.1,
    xmax = -1.1,
    ymin = 14.65,
    ymax = 16.65
  )

longitudinal_state_plot2 <- longitudinal_state_plot +
  annotation_custom(
    ggplotGrob(longitudinal_state_plot_inset),
    xmin = 1.1,
    xmax = 5.8,
    ymin = 12.65,
    ymax = 14.65
  )

leg <- get_legend(longitudinal_state_plot)

leg_plot <- as_ggplot(leg)

design <- "
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  12
  13
  13
"

longitudinal_plot <- longitudinal_rebel_plot2 + longitudinal_state_plot2 + leg_plot +
  plot_layout(design = design) &
  theme(legend.position = "none",
        plot.margin = margin(3, 3, 1.75, 2.5, "pt"))

ggsave(
  plot = longitudinal_plot,
  "figures/coefficient plots/longitudinal.pdf",
  device = cairo_pdf,
  width = 10,
  height = 6.5
)

# Longitudinal Analysis (Simple) -----------------------------------------------

## Data ------------------------------------------------------------------------

longitudinal_rebel <- read_dta("./results/longitudinal_rebel.dta") %>%
  filter(parm %in% "secessionistconflict_duration") %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict_duration" ~ "Secessionist Conflict ×\nConflict Duration"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  ) %>%
  add_row(new_parm = "Constitutive Terms &\nOther Covariates",
          sig = 1,
          model_name = "ZIOP (Outcome)")

longitudinal_state <- read_dta("./results/longitudinal_state.dta") %>%
  filter(parm %in% "secessionistconflict_duration") %>%
  mutate(
    new_parm = case_when(
      parm == "secessionistconflict_duration" ~ "Secessionist Conflict ×\nConflict Duration"
    ),
    model_name = case_when(
      eq %in% c("rebelsv_max", "statesv_max") ~ "ZIOP (Outcome)",
      eq == "inflate" ~ "ZIOP (Inflation)"
    ),
    sig = ifelse(p < 0.1, 1, 0.2)
  ) %>%
  add_row(new_parm = "Constitutive Terms &\nOther Covariates",
          sig = 1,
          model_name = "ZIOP (Outcome)")

longitudinal_rebel$new_parm <- 
  factor(longitudinal_rebel$new_parm, levels = rev(unique(longitudinal_rebel$new_parm)))

longitudinal_state$new_parm <- 
  factor(longitudinal_state$new_parm, levels = rev(unique(longitudinal_state$new_parm)))

## Plot ------------------------------------------------------------------------

space_between_bars <- 0.1

longitudinal_rebel$type_y_adj =
  scale(as.numeric(as.factor(longitudinal_rebel$model_name))) * space_between_bars

longitudinal_rebel$y =
  as.numeric(as.numeric(as.factor(longitudinal_rebel$new_parm)) + longitudinal_rebel$type_y_adj)

longitudinal_rebel_plot <- longitudinal_rebel %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  annotate(
    "text",
    x = 0,
    y = 1,
    label = paste("🗸"),
    color = "black",
    alpha = 1,
    size = 8
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)))
  ) +
  facet_wrap(paste0("Rebel SV\nObservations = 626") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:2,
    labels = rev(unique(longitudinal_rebel$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off")

space_between_bars <- 0.1

longitudinal_state$type_y_adj =
  scale(as.numeric(as.factor(longitudinal_state$model_name))) * space_between_bars

longitudinal_state$y =
  as.numeric(as.numeric(as.factor(longitudinal_state$new_parm)) + longitudinal_state$type_y_adj)

longitudinal_state_plot <- longitudinal_state %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_segment(
    aes(
      x = min90,
      xend = max90,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  annotate(
    "text",
    x = 0,
    y = 1,
    label = paste("🗸"),
    color = "black",
    alpha = 1,
    size = 8
  ) +
  guides(
    color = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    shape = guide_legend(
      reverse = TRUE,
      order = 1,
      byrow = TRUE
    ),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), byrow = TRUE)
  ) +
  facet_wrap(paste0("State SV\nObservations = 620") ~ ., scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format(
                       format = function(x)
                         as.character(x)
                     )) +
  scale_y_continuous(
    breaks = 1:2,
    labels = rev(unique(longitudinal_state$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.1)"), "No (p \u2265 0.1)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23, 22)) +
  scale_color_manual(values = c("blue", "red")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    axis.text = element_text(color = "black"),
    axis.text.y = element_blank(),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    legend.direction = "vertical",
    plot.margin = margin(0, 0, 5, 0, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(-5, 125.7, 0, -115.7, "pt"),
    legend.spacing.x = unit(5, "pt"),
    legend.spacing.y = unit(-5, "pt"),
    legend.key.width = unit(25, "pt"),
    legend.key.size = unit(12.5, "pt"),
    legend.title = element_text(size = 9, hjust = 0.5),
    legend.text = element_text(size = 8),
    legend.title.position = "left",
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"), lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off")

longitudinal_simple_plot <- longitudinal_rebel_plot + longitudinal_state_plot &
  theme(plot.margin = margin(3, 3, 1.75, 2.5, "pt"))

longitudinal_simple_plot

ggsave(
  plot = longitudinal_simple_plot,
  "figures/coefficient plots/longitudinal_simple.svg",
  width = 8,
  height = 2.5
) # Save as .pdf using Inkscape: File > Save As... > PDF
